TGCN template

1 Load data

# Load data
data <- as.data.frame(t(readRDS(paste0(data_path, "ROSMAP_RINPMIAGESEX_resids.rds"))))
covs <- readRDS(paste0(data_path, "ROSMAP_RINPMIAGESEX_covs.rds"))
stopifnot(identical(covs$mrna_id, colnames(data)))

covs <- covs[, c("ceradsc", "braaksc", "cogdx", "neuroStatus")]
covs$ceradsc <- as.factor(covs$ceradsc)
covs$braaksc <- as.factor(covs$braaksc)
covs$cogdx <- as.factor(covs$cogdx)
covs$neuroStatus <- as.factor(covs$neuroStatus)

table(covs$neuroStatus)

## 
##   0   1 
## 168 174

protein <- "APP"
input <- data[-match(protein, rownames(data)), ]
toPredict <- as.vector(unlist(data[match(protein, rownames(data)), ]))

2 Cutoff selection

r <- testAllCutoffs(exprData=input,
                    target=toPredict,
                    covs=covs,
                    train.split=0.7,
                    nfolds=5,
                    t=10,
                    path=paste0(tgcn_path),
                    targetName="APP",
                    tissueName="ROSMAP",
                    seed=1234,
                    cutoffs=10:1,
                    n=100, 
                    m=10, 
                    s=10, 
                    minCor=0.3,
                    maxTol=3,
                    save=T,
                    overwrite=F)

r <- readRDS(paste0(tgcn_path, "/Net/APP_ROSMAP_TGCNs_all_cutoffs.rds"))

grid.arrange(r$selectCutoff$nHubs, r$selectCutoff$stats + theme(legend.position="bottom"))

p <- lapply(r$nets, function(cutoff) cutoff$GOenrich$plotStats) 

ggarrange(p$c8 + theme(text=element_text(size=10)) + scale_y_continuous(limits=c(0,100)), 
          p$c9 + theme(text=element_text(size=10)) + scale_y_continuous(limits=c(0,100)), 
          p$c10 + theme(text=element_text(size=10)) + scale_y_continuous(limits=c(0,100)), 
          ncol=3, nrow=1, common.legend=T, legend="bottom")

3 TGCN for each cutoff

3.1 Cutoff 8

3.1.1 Modules size selection

r$nets$c8$net$moduleSizeSelectionPlot

3.1.2 Modules correlation

r$nets$c8$net$plotCorr

3.1.3 Modules composition

DT::datatable(r$nets$c8$net$modules)

3.1.4 Cross tab plot

3.1.5 GO enrichment stats

grid.arrange(r$nets$c8$GOenrich$plotStats, r$nets$c8$GOenrich$plotNterms, nrow=2)

DT::datatable(r$nets$c8$GOenrich$terms)

3.1.6 Reduced GO terms at TGCN level

r$nets$c8$GOenrich$plotReduced

reduced_terms <- r$nets$c8$GOenrich$reducedTerms
treemapPlot(reduced_terms)

wordcloudPlot(reduced_terms, min.freq=1, colors="black")

3.1.7 Cell-type enrichment

r$nets$c8$CTenrich$plot

r$nets$c8$CTenrich$sigTests

3.1.8 Module-trait corr

r$nets$c8$moduleTraitCorr$plot_pval

3.2 Cutoff 9

3.2.1 Modules size selection

r$nets$c9$net$moduleSizeSelectionPlot

3.2.2 Modules correlation

r$nets$c9$net$plotCorr

3.2.3 Modules composition

DT::datatable(r$nets$c9$net$modules)

3.2.4 Cross tab plot

3.2.5 GO enrichment stats

grid.arrange(r$nets$c9$GOenrich$plotStats, r$nets$c9$GOenrich$plotNterms, nrow=2)

DT::datatable(r$nets$c9$GOenrich$terms)

3.2.6 Reduced GO terms at TGCN level

r$nets$c9$GOenrich$plotReduced

reduced_terms <- r$nets$c9$GOenrich$reducedTerms
treemapPlot(reduced_terms)

wordcloudPlot(reduced_terms, min.freq=1, colors="black")

3.2.7 Cell-type enrichment

r$nets$c9$CTenrich$plot

r$nets$c9$CTenrich$sigTests

3.2.8 Module-trait corr

r$nets$c9$moduleTraitCorr$plot_pval

3.3 Cutoff 10

3.3.1 Modules size selection

r$nets$c10$net$moduleSizeSelectionPlot

3.3.2 Modules correlation

r$nets$c10$net$plotCorr

3.3.3 Modules composition

DT::datatable(r$nets$c10$net$modules)

3.3.4 Cross tab plot

3.3.5 GO enrichment stats

grid.arrange(r$nets$c10$GOenrich$plotStats, r$nets$c10$GOenrich$plotNterms, nrow=2)

DT::datatable(r$nets$c10$GOenrich$terms)

` ### Reduced GO terms at TGCN level

r$nets$c10$GOenrich$plotReduced

reduced_terms <- r$nets$c10$GOenrich$reducedTerms
treemapPlot(reduced_terms)

wordcloudPlot(reduced_terms, min.freq=1, colors="black")

3.3.6 Cell-type enrichment

r$nets$c10$CTenrich$plot

## NULL

r$nets$c10$CTenrich$sigTests

## NULL

3.3.7 Module-trait corr

r$nets$c10$moduleTraitCorr$plot_pval

4 Hub gene module enrichment

go_results <- r$nets$c8$GOenrich$terms
plots <- getReducedTermsPlots(go_results, module=T)
saveRDS(plots, paste0(tgcn_path, "plots_APP.rds"))

go_results <- r$nets$c8$GOenrich$terms
sort(table(go_results$query), decreasing=T)

## 
##     RPS10 HNRNPA2B1    ATP1B1    PPP2CA     RPL36      TPPP    FBXO11     SCN8A 
##       215        85        79        66        63        55        53        53 
##     NFASC       NCL     GNAI1     ABCD3     PTAR1     YWHAE     DEAF1   SERINC3 
##        52        44        22        16        16        14        12        11 
##      FRZB     FBXO9    CEP250    DGCR6L     SCAF4      hubs   PRKAR1A     ZFPM1 
##        11         6         5         5         3         3         1         1 
##     ASTE1     PLPP5 
##         1         1

plots <- readRDS(paste0(tgcn_path, "plots_APP.rds"))

4.1 ATP1B1

plots$ATP1B1$scatterPlot

4.2 DEAF1

plots$DEAF1$scatterPlot

4.3 FBXO11

plots$FBXO11$scatterPlot

4.4 HNRNPA2B1

plots$HNRNPA2B1$scatterPlot

4.5 NCL

plots$NCL$scatterPlot

4.6 NFASC

plots$NFASC$scatterPlot

4.7 PPP2CA

plots$PPP2CA$scatterPlot

4.8 RPL36

plots$RPL36$scatterPlot

4.9 RPS10

plots$RPS10$scatterPlot

4.10 TPPP

plots$TPPP$scatterPlot